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ABSTRACT 

Two  methods  of  obtaining  sensitivity  data  were  simulated  on  an 
electronic  computer  for  the  purpose  of  comparing  the  accuracy  of  the 
estimates  of  the  parameters  of  an  underlying  cumulative  normal  response 
function.   The  first  method  simulated  the  standard  Bruceton  procedure 
while  the  second  used  a  modified  binary  search  routine  with  a  portion  of 
the  sample  in  order  to  obtain  maximum  likelihood  estimates  of  the  input 
parameters  for  use  in  a  follow-on  Bruceton  test. 

The  results  showed  both  methods  to  be  effective  in  estimating  the 
mean  but  with  slightly  more  variability  in  the  estimates  obtained  by  the 
second  procedure.   Both  methods  underestimated  the  standard  deviation  - 
again  with  more  variability  in  the  estimates  obtained  by  the  second 
procedure.   When  the  prior  parameter  estimates  were  unknown  and  the 
applicable  stimulus  level  bounded,  the  second  method  yielded  estimates 
favorably  comparable  to  those  expected  from  the  Bruceton  procedure  with 
suitable  prior  input  estimates. 
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I.   INTRODUCTION 

Frequently  a  statistician  is  faced  with  the  problem  of  determining 
the  level  of  a  stimulus  which  critically  affects  the  performance  of  a 
device.   The  nature  of  the  testing  to  be  discussed  is  such  that  once 


some  positive  level  of  the  stimulus  is  applied  to  the  device  either  a 
response  or  a  non-response  can  be  immediately  observed  and,  in  either 
case j  the  device  is  altered  so  that  a  bonafide  result  cannot  be  obtained 
from  a  second  test.   Tests  of  this  type  are  known  as  sensitivity  tests. 

One  of  the  many  problems  besetting  those  involved  in  explosives 
research  is  that  of  providing  measures  and  specifying  rules  to  provide 
for  the  safe  handling  and  transportation  of  explosives.   Many  different 
types  of  sensitivity  testing  apparatus  have  been  developed  for  laboratory 
use,  the  most  common  being  those  that  subject  some  quantity  of  explosive 
to  the  impact  load  of  a  falling  drop-weight  from  some  controllable  height. 
At  least  as  late  as  October  1965  there  remained  two  important  physical 
problems  to  be  solved;  namely,  that  of  establishing  a  measure  of  stimulus 
not  highly  apparatus -dependent  and  then  that  of  translation  of  these 
results  to  safe  handling  rules  [1].   These  problems  are  not  addressed  in 
this  paper  but  should  be  kept  in  mind  when  considering  the  overall  problem, 

In  the  early  1940's,  a  technique  for  obtaining  sensitivity  data  was 
developed  and  used  in  explosives  research  at  the  Explosives  Research 
Laboratory,  Bruceton,  Pennsylvania  which  has  come  to  be  called  synony- 
mously, the  Bruceton,  Staircase,  or  "Up  and  Down"  Method. 

The  aim  of  this  method  of  testing  is  to  increase  the  accuracy  with 
which  certain  critical  values  of  the  stimulus  may  be  estimated,  notably 
the  median  (or  mean)  and  standard  deviation.   The  accuracy  of  the  method 


depends  in  part  on  the  stimulus  level  at  which  the  first  item  is  tested 
and  the  interval  spacing  for  subsequent  levels  of  testing  [2]. 

When  the  stimulus  levels  mentioned  above  cannot  be  determined  prior 
to  testing  or  when  little  confidence  is  placed  on  the  available  estimates, 
a  preliminary  (or  search)  phase  of  testing  may  be  desirable  to  obtain 
maximum  likelihood  estimates  prior  to  employing  the  Bruceton  Method  with 
the  remainder  of  the  sample.   A  procedure  to  do  this  is  offered  as  an 
alternative  method. 

The  comparative  accuracies  of  the  two  techniques  were  examined 
through  the  use  of  simulation  conducted  on  a  high-speed  electronic 
computer.   All  parameters  and  estimates  considered  as  inputs  to  the 
simulation  were  kept  within  ranges  for  which  the  Bruceton  Method  is 
considered  to  yield  accurate  results  [2]. 


II.   THE  MODEL 

Let  x  be  an  applied  stimulus  level  (xe;[o,<»))  and  y  =  y(x)  be  the 
associated  response  (ye<o,l>  where  "o"  denotes  no  response  and  "1" 
denotes  response).   At  any  given  stimulus  level  consider  y  to  be  the 
realization  of  a  Bernoulli  random  variable,  Y,  with  response 
probability 

p(x)  =  Prob  (Y  =  l|x) 
The  function  p(x)  is  called  the  response  function  and  is  further 
specified  as 

p(x)  =  0  xe[o,a] 

0  <  p(x)  <  1        xe(a,b) 
and  p(x)  =  1  xe[b,=°) 

The  intervals  [o,a],  (a,b)  and  [b,°°)  are  called  the  zero-response 
region,  the  mixed-response  region,  and  the  one-response  region 
respectively.   It  is  assumed  that  p(x)  is  a  monotonely  increasing 
function  for  stimulus  values  in  the  mixed-response  region.   Thus, 
p(x)  can  be  considered  as  the  cumulative  distribution  function  for  a 
random  variable  X  such  that 

p(x)  =  Prob  (X  <  x).   [3] 
In  this  context  the  random  variable  X  can  be  interpreted  as  a  thres- 
hold stimulus  level,  thus 

Prob  (Y  =  l|x)  =  Prob  (X  <  x)  =  p(x) 

and  Prob  (Y  =  o|x)  =  Prob  (X  >  x)  =  1  -  p(x).   [3] 

2 

It  is  assumed  the  X  is  distributed  Normal  (\a,o   );  that  is 

p(x)  =  Cp(x|(i,0  ) 
where  cp(x^,a  )  represents  the  cumulative  normal  distribution  with  mean 


2  .    , 

|i  and   variance   a   .       In   particular 


Prob    (x  <  \a)   =  p(\d)   =0.5.      [3] 
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III.   TESTING  METHODS 

A.   BRUCETON  METHOD 
1 .   Description 

Based  on  intuition  or  past  experiments,  the  experimenter  selects 
a  priori  estimates  of  p.  and  a.   Call  these  estimates  |j  and  o_  and  let 


d  =  or 


The  experimenter  tests  the  first  item  at  or  near  n  .   If  there 


is  a  response  the  second  item  is  tested  at  a  level  d  units  below  \i    , 
otherwise  the  second  item  is  tested  at  a  level  d  units  above  \i    .   In 
the  same  manner,  each  of  the  remaining  items  is  tested  at  a  level  d 
units  above  or  below  the  previous  test  level  according  as  there  was  not 
or  there  was  a  response  observed  for  the  previous  test.   Thus  the 
sample  is  concentrated  about  the  mean  and  one  would  expect  nearly  equal 
numbers  of  responses  and  non-responses.   In  fact,  the  number  of  non- 
responses  at  any  level  will  not  differ  by  more  than  one  from  the 
number  of  responses  at  the  next  higher  level  [2]. 

Let  N  denote  the  total  number  of  observations  of  the  less 
frequent  event  and  nft,n..  ,n~  ,  •  •  tl  denote  the  frequencies  of  this  event 
at  each  level  where  n_  corresponds  to  the  lowest  level  and  n,  the 
highest  level  at  which  the  less  frequent  event  occurs. 

The  final  estimates  of  \A   and  a  are  based  on  the  first  two 

moments  of  the  stimulus  levels.   Since  the  intervals  are  equally 

spaced,  these  moments  can  be  computed  in  terms  of  the  sums 

A  =  £  i  n. 
1 

l 

2 

and  B  =  S  i  n. . 

i     X 
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Let  [A   be  the  estimate  of  g.  by  this  method.   Then 


&-«•  +  d(|±D 


where  x'  represents  the  lowest  level  at  which  the  less  frequent  event 
occurs  [2].  The  plus  sign  is  used  when  the  analysis  is  based  on  non- 
responses,  and  the  minus  sign  when  it  is  based  on  responses  [2]. 

2   2 

If  (NB-A  )/N  >  .3  the  sample  standard  deviation  is 


2 

s  =  1.620  d  (~p_  +  .029^) 


m 


2 

zl 

o 

N 

Otherwise,  a  more  elaborate  calculation  must  be  employed  and  is 

described  in  Ref.  2. 

To  obtain  confidence  intervals,  estimates  of  the  standard 

deviations  of  the  sample  mean  and  sample  standard  deviation,  say  s 

and  s   respectively,  are  given  by 
s 

_  Gs 

and 

Hs 
s   =  — 

s  J* 

where  the  factors  G  and  H  are  dependent  on  the  ratio  —  and  the 

r  s 

position  of  the  mean  relative  to  the  testing  levels.   Plots  of  these 
factors  are  available  in  Ref.  1. 
2 .   Discussion 

Only  rarely  is  the  threshold  stimulus  Z  normally  distributed. 
It  is  usually  the  case  that  some  scale  transformation  of  Z,  say  X,  is 
made  so  that  X  is  normally  distributed  in  the  vicinity  of  the  mean. 
This  transformation  is  done  prior  to  testing  to  determine  \±     and  a  . 
Only  after  all  analysis  is  completed  are  the  values  scaled  back  to  the 
original  stimulus  measure  [2]. 

12 


The  size  of  the  sample  is  critical  to  the  accuracy  of  the 
estimation.   Note  that  at  most  only  half  of  the  sample  is  used  in  the 
analysis  so  that,  for  example,  if  thirty  items  are  tested  the  maximum 
possible  value  of  N  is  fifteen.   The  analysis  is  based  on  large  sample 
theory  which  in  the  case  mentioned  would  be  applied  to  a  sample  of  size 
fifteen  [2]  [4]. 

Unless  normality  of  the  variate  is  assured  this  method  does 
not  yield  accurate  results  for  the  small  end    large  percentage  points. 
This  is  unfortunate  since  in  most  applications  one  would  be  more 
interested  in  a  small  percentage  point  as  a  measure  of  safety  and  a 
large  percentage  point  as  a  measure  of  reliability.   At  any  rate,  an 
estimate  of  a  percentage  point  j  is 

A 

j  =  n  +  ks 
where  k  is  chosen  from  tables  of  the  standard  normal  deviate  to  give 
the  desired  percentage  [2].   One  could  then  conduct  tests  in  the 
vicinity  of  this  value  to  refine  the  estimate. 

B.   BRUCETON  METHOD  PRECEDED  BY  SEARCH 
1.   Description 

In  the  event  that  a  priori  estimates  of  [i   and  a  are  not 
available  some  economic  method  of  attaining  these  estimates  is  desired. 
A  method  proposed  and  described  below  is  a  modified  binary  search 
technique. 

Again,  the  assumption  is  that  the  threshold  stimulus  (or 
some  transformation  of  it)  is  normally  distributed  and  p(x)  can  be 
represented  by  a  cumulative  normal  distribution. 

As  noted  from  the  model 

Prob  (Y  =  0|x  <  a)  =  1 
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and 

Prob  (Y  =  l|x  >  b)  =  1. 
The  first  step  in  the  procedure,  then,  is  to  select  values  for  a  and  b. 
(In  the  case  of  complete  uncertainty  these  could  be  the  limiting  values 
of  the  testing  apparatus)  and  commence  the  binary  search  starting  at 

x  =  (a  +  b)/2. 
If  p(x)  were  a  step  function,  repetition  of  this  method  would  locate 
the  step  in  an  interval  of  any  desired  length.   In  general,  however, 
the  mixed-response  region  has  non-zero  width  and  a  non-response  would 
merely  indicate  that  the  applied  stimulus  is  in  the  mixed  response 
region  or  below  while  a  response  would  indicate  that  it  was  in  the 
mixed  response  region  or  above. 

If  a  test  at  x1  yields  a  response  and  a  test  at  x  yields  a 
non-response  while  x..  <  x„  it  is  certain  that  both  x-  and  x?  are  in  the 
mixed  response  region.   This  condition  is  called  a  response  inversion 
and  is  the  basic  indicator  for  the  modified  binary  search  technique. 
The  description  of  the  procedure  is  best  followed  by  referring  to 
Figures  1  through  4. 

Sequence  S*  is  a  cyclic  one  indicating  that  a  reduction  in  step 
size  should  be  taken.   Test  levels  are  selected  attempting  to  reproduce 
this  sequence.   Failure  to  do  this  results  in  the  basir  inversion 
sequence  S  .   Tests  are  then  made  at  the  end  of  this  sequence  to  result 
in  one  of  three  terminal  situations  S1 ,  S~  ,  or  S_.   In  the  event  the 
mixed  response  region  is  relatively  narrow  and  near  a  or  b ,  several 
binary  reductions  may  be  necessary  to  reproduce  S*  or  one  of  the 
terminal  situations.   These  circumstances  are  represented  by  S   and 
S„  [31. 
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STARTING  SEQUENCE 
Figure  1 
15 


yg^» 


v«^ 


r^s 


s, 


S      SEQUENCE 
Figure   2 


16 


p      O        I       » 
£      **.    *j    if/ 


*£ 


*x+*3 


AK\ 


U& 


WO 


y-,-2^-^ 


Y<?*s 


V-,  =2K,-yc 


Aj*> 


yg^. 


S2. 


*l 


S*  SEQUENCE 
Figure   3 
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yes 


wo 


JVL 


K>d> 


S      SEQUENCE 
Figure   4 
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Maximum  likelihood  estimates  of  \i   and  a  are  available  for 
sequences  S1 ,  S_,  and  S_  and  developed  as  described  below  [3]. 


2.   Discussion 


It  is  assumed  that  all  trials  are  independent.   Thus  the 
probability  of  the  sequence  S,  is 
Prob  (sp  -  Prob  (Y^O,  Y2=0,  Y3=l,  Y4=0,  Y5=l,  Y&=1  |x][  ^  ,X3  ,X4,X5  ,Xfi) 


where 


=  it  Prob  (Y.  =  y.  x.) 
.  ,        i    i   i 
1=1 


Prob  (Y.  =  y.  x.)  =  cp(x.)  if  y.  =  1 


and 


=  1  -  cpCx^  if  y  =  0 


cp(x.)  =  Prob  (X.  -:  x.)  =  f  — - —  e  2    ^       dx. 

1_   L     J   72TTa 

-00 

Maximum  likelihood  estimates  for  \±   and  a  can  then  be  established 
using  standard  normal  tables  for  each  of  the  terminal  situations. 
These  estimates  are  indicated  on  Figure  5. 
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RESPONSE   SEQUENCES   AND   PARAMETER  ESTIMATES 

Figure   5 
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IV.   SIMULATION 


A.   DESCRIPTION 

All  simulated  experiments  were  conducted  on  an  IBM  360/67  computer 
using  the  FORTRAN  IV  programming  language.   The  basic  program  is 
attached.   The  response  function  p(x)  used  was  cumulative  normal  with 
[X   =  30  and  a  =  3 . 

The  sample  size  was  kept  at  seventy  for  each  experiment  to  provide 
some  assurance  that  the  analytical  sample  would  be  suitable  for  large 
sample  analysis. 

The  basic  test  procedure  was  to  draw  a  random  number  on  the  unit 
interval  and  compare  this  to  F(x),  a  function  of  a  standard  normal 
variate  specified  as 


and 


F(x)  -  |  [l  -  erf  Q0\ 


F(x)  -  I  [l  +  erf  (jp[ 


if  x  <  0 


if  x  >  0 


where 


erf (v)  =~   e 


2 


dt 


0 


(The  function  subprogram  erf  is  an  IBM-supplied  subprogram.)   If  the 
random  number  was  less  than  or  equal  to  F(x.)  then  a  response  was 
counted  for  the  i    level;  otherwise  a  non-response  was  counted. 

Six  different  cases  were  tested  using  the  straight  Bruceton  pro- 
cedure (METHOD  1)  with  two  different  input  estimates  of  \i   and  three 
different  input  estimates  of  a.   Case  1  considered  exact  estimates; 
i.e.,  (j  =  [l   and  <j     =  a.   Case  2  considered  u  =  ^i-6  and  a,  =  a- 
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Cases  3  and  4  considered  |j  =  |j  and  aT  =  a/2,  2<j  respectively  while 
Cases  5  and  6  repeated  Cases  3  and  4  except  |j.  =  (j-6.   For  each  of  the 
six  cases  1000  experiments  were  conducted  each  utilizing  a  different 
sequence  of  random  numbers . 

The  search  procedure  (METHOD  2)  was  then  incorporated  into  each 
of  the  above  six  cases  using  the  a  prior  estimates,  \i     and  oT>  to 
determine  estimates  for  stimulus  levels  a  and  b  and  thereby  the  size 
of  the  binary  reduction  as  indicated  in  Figure  1.   The  program  then 
followed  the  flow  shown  in  Figures  1  through  4  until  either  a  terminal 
sequence  was  reached  or  the  search  was  arbitrarily  terminated  as  dis- 
cussed in  subparagraph  C  below.   The  Bruceton  procedure  was  then  used 
until  the  sample  was  exhausted. 

The  final  case,  Case  7,  indicated  complete  lack  of  knowledge  of  (j. 
and  a  but  considered  the  upper  and  lower  stimulus  level  limits  of  the 
test  apparatus  to  be  100  and  0  respectively. 

B.   MEASURES  OF  EFFECTIVENESS 

At  the  completion  of  all  experiments  for  each  case,  several 
measures  were  obtained  for  comparison.   First,  average  values  of  the 
parameters  were  determined  to  be 

A       A   , 

\X   -  S  l^/N 
i 

and 

a     A  , 
a  =  E  aVN 

i 

A      A  th 

where  \j.     and  o  are  the  a  posteriori  estimates  of  \jl   and  q  for  the  i 

experiment  and  n,  the  number  of  experiments  used.   Next,  as  measures  of 

variability 
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and 


2       a      2 

s  a  -  S  ((j  -  |_i)  /n-1 

^   i 


s  a  =  S  (a.  -  a)/n-l 
a   .   i 


were  calculated.   In  addition,  the  program  listed  the  maximum  and 
minimum  estimates  of  both  [j  and  a. 

C.   DISCUSSION 

In  Chapter  III  it  was  noted  that  sequences  S*,  S  ,  and  S   are 

U         Li 

cyclic.   In  order  to  simplify  the  program  it  was  necessary  to 
artificially  terminate  these  situations  at  some  point  and  calculate 
the  input  values  for  the  Bruceton  test.   The  estimate  of  \j.   used  was 

where  x1  and  x_  are  adjacent  testing  levels  and  x9  >  x,  with  y  =  0 
and  y_  =  1.   The  estimate  of  a  used  was 


y2 


as  -  (x2  -  xx)/2 


for  Cases  1  through  6  and 

°s  =  (x2  •  xl)/6 
for  Case  7.   The  former  estimate  of  a  was  chosen  arbitrarily  while  the 
latter  estimate  was  based  on  the  estimate  of  the  mixed  response  region 
being  6a.   While  the  number  of  terminations  of  this  type  was  insignifi- 
cant for  the  first  six  search  cases,  in  the  final  case  over  600  experi- 
ments were  thus  terminated  requiring  the  program  to  be  expanded  to 
permit  more  recycling.   The  point  is  that  the  artificial  termination 
does  not  represent  the  search  procedure.   This  problem  would  not  arise 
in  field  experimentation  until  either  the  sample  was  exhausted  or  the 
step  size  reduction  of  stimulus  level  indicated  was  too  narrow  to  be 

measured  or  controlled  by  the  test  apparatus. 
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Also  in  the  interest  of  program  simplification  those  experiments 
for  which 

F__L_a!  <     3 
2     -  ,J 
N 

were  not  used  for  analyses.   This  limitation  invalidated  the  measures 

of  effectiveness  for  the  Bruceton  cases  where  a  =  2<j. 

D.   RESULTS 

The  results  of  the  simulation  are  listed  in  Table  I.   It  is 
questionable  that  the  measures  listed  under  Method  1  are  valid  for 
Cases  4  and  6  in  that  only  .381  and  .393  of  the  possible  experiments 
were  used.   These  two  cases  and  Case  4  under  Method  2  (where  .661  of 
the  possible  experiments  were  used)  are  the  only  ones  for  which 

O   >    O. 

In  general  the  extreme  estimates  are  more  widely  separated  and  the 
variability  of  a  is  greater  in  Method  2. 

Estimates  of  \j.   range  from  27.8823  to  31.7647  for  Method  1  and 
27.937  to  31.91  for  Method  2. 

Estimates  of  a  range  from  .8741  to  6.5027  for  Method  1  and  .3498 
to  9.8328  for  Method  2. 

A 

The  lowest  average  \l,   29.9113,  was  obtained  under  Method  1,  Case  5, 
while  the  highest  average  |j,  30.1175,  was  obtained  under  Method  2, 
Case  3. 

The  lowest  average  a,  2.3748,  was  obtained  under  Method  2,  Case  5, 
while  the  highest  average  a,  2.9474,  was  obtained  under  Method  1, 
Case  5.   (Case  6  is  not  counted  under  Method  1  nor  is  Case  4  under 
both  methods . ) 
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TABLE   OF   EXPERIMENTAL  RESULTS 


METHOD    1 

METHOD   2 

A 

M 

A 

a 

A 

M 

A 

a 

CASE 

1 

Mj  =  30 

AVE 

30.0067 

2.8320 

30.0117 

2.8609 

ax  =  3 

MAX 

31.7647 

5.7904 

31.7813 

5.9343 

a  =   18 

MIN 

28.5000 

1.6089 

28.2187 

1.1241 

b  =  42 

VAR 

.2523 

.4128 

.2514 

.5831 

CASE 

2 

Mj  =  24 

AVE 

29.9641 

2.9040 

30.0317 

2.8819 

ax  =  3 

MAX 

31.6765 

5.8249 

31.6875 

9.1369 

a  =   12 

MIN 

28.3235 

1.6250 

28.1976 

.9512 

b  =  36 

VAR 

.2656 

.4225 

.2666 

.6336 

CASE 

3 

Mx  =  30 

AVE 

30.0295 

2.7216 

30.1175 

2.8615 

a,  =  1.5 

MAX 

31.6071 

6.1197 

31.9100 

7.7997 

a  =  24 

MIN 

28.4118 

.8741 

28.5950 

.8697 

b   =  36 

VAR 

.2046 

.7409 

.2693 

.9081 

CASE 

4 

Hj  =  30 

AVE 

29.9683 

3.5424 

29.9750 

3.0721 

a,  =  6 

MAX 

31.4571 

6.0266 

31.6875 

6.3569 

a  =  6 

MIN 

28.0286 

-- 

27.9370 

1.6170 

b  =  54 

VAR 

.2574 

.4639 

.2619 

.4522 

CASE 

5 

Hj  =  24 

AVE 

29.9113 

2.9474 

29.9363 

2.3748 

Qj  =   1.5 

MAX 

31.4773 

5.9257 

31.4063 

7.9507 

a  =   18 

MIN 

28.2353 

.9452 

28.1961 

.3498 

b  =  30 

VAR 

.2220 

.8748 

.2184 

1.4889 

CASE 

6 

Mj   =  24 

AVE 

29.9493 

3.5438 

30.0247 

2.8398 

a  =  6 

MAX 

31.4118 

6.5027 

31.5756 

6.4201 

a  =  0 

MIN 

27.8823 

-- 

28.2552 

1.1252 

b  =  48 

VAR 

.2639 

.4785 

.2490 

.6300 

CASE 

7 

-- 

AVE 

-- 

-- 

30.0123 

2.7280 

-- 

MAX 

-- 

-- 

31.8229 

9.8328 

a  =  0 

MIN 

-- 

-- 

27.9541 

.5082 

b  =   100 

VAR 

-- 

-- 

.2628 

2.1041 

TABLE    I 
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V.   CONCLUSIONS  AND  RECOMMENDATIONS 

A.  CONCLUSIONS 

1.  Estimation  of  the  Mean 

Both  methods  estimate  the  mean  effectively. 

2 .  Estimation  of  the  Standard  Deviation 

Both  methods  tend  to  under-estimate  the  standard  deviation  with 
no  predictable  bias  and  are  therefore  unsuitable  for  use  in  safety  or 
reliability  statements.   This  conclusion  agrees  with  the  findings  of 
Hampton  [4]  as  it  pertains  to  the  Bruceton  Method. 

3 .  Extension  of  the  Search  Phase  for  the  Starting  Sequence 
Termination  of  the  search  phase  with  sequence  S,  in  the  starting 

sequence  (see  Figure  1)  may  yield  estimates  of  a  greater  than  twice 
the  actual  value.   To  avoid  this  it  is  advisable  to  extend  the  search 
phase  as  described  in  Ref.  3. 

4.  Use  of  Search  Technique 

The  search  procedure  should  be  used  in  those  cases  where  there 
is  not  independent  evidence  that  the  estimate  of  a  is  within  the  range 
for  which  the  Bruceton  Method  is  recommended  (i.e.,  a/2  <  aT  <  2a). 

B.  RECOMMENDATIONS 

Further  testing  of  Method  2  is  recommended  under  the  circumstances 
listed  below. 

1.   Reduction  of  Sample  Size 

It  would  be  of  interest  to  reduce  the  sample  size  to  the  point 
where  the  effective  sample  is  small,  say  15,  and  compare  the  Bruceton 
procedure  with  the  search  procedure  using  the  entire  sample  for  the 
search. 
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2 .   Random  Selection  of  Response  Function  Parameters 

A  more  valid  test  of  both  methods  would  be  achieved  by 
randomally  selecting  values  of  \j   and  a   over  some  range  and  using 
the  on-line  computer  facility  to  conduct  the  simulation. 
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COMPUTER  PROGRAM 

THIS  PROGRAM  SIMULATES  SENSITIVITY  TESTING  BY  BOTH  THE  RRUC- 
FTQNI  METHOD  (WHEN  IANY  =  0)  AND  THF  BRUCETON  METHOD  PRFCEDED 
BY  THE  MODIFTHI  BINARY  SEARCH  (WHT-N  IANY=1).THE  UNDERLYING 
RESPONSE  FUNTION  IS  CUMULATIVE  NORMAL  C*0,  "*)  .THP  IN^UT  EST- 
IMATES OF  THE  MEAN  AND  TH=  STANDARD  DEVIATION  ARE  CALLED 
EXMij  AND  EXSIG  RESPECTIVELY. 

THE  PRINCIPLE  VARIABLE  NAMES  ARF  AS  FOLLOWS... 

AACT  IS  THP  STIMULUS  VALUE  AT  THE  UPPER  LIMIT  OF  THE 

MIXED  RESPONSE  REGION. 
BACT  IS  THC  STIMULUS  VALUE  AT  THE  LOWER  LIMTT  OF  THE 

MIXED  RESPONSE  REGION. 
A  AND  B  ARE  ESTIMATES  OP  AACT  AND  RACT  RESPCCT I VFLY. 
X(J)  IS  THE  STIMULUS  LEVEL  OF  THE  JTH.  STIMULUS. 
IXO(J)  IS  THE  CUMULATIVE  COUNT  OF  NON-RESPONSES  AT  X(J). 
IXX(J)  TS  THF  CJJMUIATIVF  COUNT  OF  RESPONSES  AT  XU) 
IS  IS  THE  SAMPLE  SI7E. 
NU  IS  THE  ENTRY  NUMBER  FOR  THE  RANHOM  NUMRER  GENERATOR, 

UN  I  F  . 
N  COUNTS  THE  NUMBER  OF  EXPERIMENTS. 

RN  IS  THE  RANDOM  NUMBER  ON  (0,1)  RETURNED  RY  UNIF. 
FOEX  IS  THE  VALUE  OF  THE  RESPONSE  FUNCTION  RETURNFO  BY 

SUBPROGRAMS  XNCDF  AND  SNCDF. 
I  SUMO  TS  THE  TOTAL  NUMRFR  OF  NON-RESPONSES  FOR  ONE  EXPER- 
IMENT. 
ISUMX  IS  THC  TOTAL  NUMBER  OF  RESPONSES  FOR  ONE  EXPERIMENT 
NT  IS  THE  MINIMUM  of  TSUMO  AND  ISUMX. 

NS(J)  IS  THE  FPEOUcNCY  HF  THE  LESS  FREQUENT  EVENT  AT  X(J) 
NG(J)  REARRANGES  NS(J)  SO  THAT  NG(1)=NS(I)  WHERC  X(I)  IS 

THE  LOWEST  STIMULUS  LEVEL  AT  WHICH  THE  LESS  FREQUENT 

EVENT  OCCURS. 
AR(J)  IS  USED  to  CALCULATE  THE  first  MOMENT ,SUMAP . 
BR(J)  IS  USED  TO  CALCULATC  THE  SFCOND  MOMENT, SUMRR . 
YPRIME  IS  THE  LOWEST  LEVEL  AT  WHICH  THE  LESS  FREQUENT 

EVENT  OCCURS 
XMUEST  IS  THE  FINAL  ESTIMATE  Oc  THF  TRUE  M-<\N,XMU. 
DEVEST  IS  THE  FINAL  ESTIMATE  OF  THE  TRUE  STANDARD  DEVIAT- 

EMUuVlS  THF  oiFFEPENCE  OF  XMUEST  AND  XMU. 

EDEVU)  IS  THE  DIFFERENCE  OF  DEVEST  AND  XSIG. 

SAMAVM  AND  SAMAVD  ARE  THE  SAMPLE  AVERAGE  FPRORS  OF  XMlJFST 

AND  DEVEST  R E SPECTI VFLY. 
SAMSQM  AND  SAMSQD  ARE  THE  AVERAGE  mc^n  SQUARE  errors  OF 

XMUEST  AND  DEVEST  RESPECTIVELY. 
NOGH  IS  THE  NUMRER  OF  EXPERIMENTS  NOT  USED  IN  THE  FINAL 

ANALYSIS 

DIMENSION  ARRAYS  AND  FORMAT 

SIMULATE  RRUCETON  FIRST  THEN  SEARCH 
IANY=0 
69    THING=0. 

INITIALIZE  INTERNAL  AND  OUTPUT  VARIABLES 

SET  EXPERIMENT  COUNTER,  SAMPLE  SIZE  COUNTcR,  ANH  NU. 

LC0UNT=1 
NU=1??71 
103   N=l 

IF(  IANY.FQ.O)  MBR  =  0 
IF(  IANY.EQ.l)  NBR  =  1 

SET  INPUT  VARIABLES 
XMU=30. 
XSIG=3. 
A  =  0. 
B=100. 
EXMU-50. 
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FXSIG=12.5 
A=  FXMU-IQ*FXSTO 
B  =  EXMU*IQ*FXSIG 
Xl=( A+B) /2. 

INC  =  0 

PROVIDE  BRANCH  TO  STANDARD  RRUCETON 

IF<NBR.EO.O)  GO  TO  33 

CONDUCT  SFARCH 


CALL  UNIF 
FOFX=XNCD 
IFCRN.LE. 
X2=( P+Xl) 
NBR=NBR+1 
CALL  UNIP 
FOFX=XNCO 
IF(«N.LE. 
X3  =  < R+X2) 
NBR=NBR+1 
CALL  UNIF 
FOFX=XNCD 
IFCRN.LE. 
X4=(B+X3) 
NBR=NBR+l 
CALL  UNIF 
FOFX=XNCD 
IFfRN.LF. 
X5=(R+XA) 
NBR=NBR+1 
CALL  UN' IP 
FOFX=XNCD 
IFCRN.LE. 
EXMU=( B+X 
FXSIG=(X5 
EX SI  0=2.* 
EXSTG=EXS 
GO  TO  700 

1313  FXMU=(X5+ 
EXSIG=<X5 
EXSIG=2.* 
EXSIG=EXS 
GO  TO  7C0 

9063  X5=(X3+X2 
NBR=NRP+1 
CALL  UNIF 
FOFX=XNCD 
IF(RN.L£. 
X6=< X3+X4 
NBR=NRR+1 
CALL  UNIF 
FOFX=XNCD 
IF(RN.LE. 
EXMU=( X6  + 
EXSIG=(X6 
EXSI0=2.* 
EXSIG=-XS 
GO  TO  700 

1316  FXMU=( X6+ 
EXSIG=( X6 
EX SI  0=2.* 
EXSIG=FXS 
GO  TO  700 

1314  X6=2.*X2- 
NBR=NBP+1 
CfiLL  UNIF 
FOFX=XNCD 
IFCRN.LE. 
XB=X5 


(RN,NU) 
F(X1 .X^U 
FO^X)  GO 
/2. 

{PN,NM) 
F(  X?,XMI) 
FO^X)  GO 
/2. 

(RNtNU) 
F( X3TXMU 
FD^X)  GO 
/2. 

(RNtNU) 
F( XA,X^U 
FOFX)  GO 
/2. 

CRN,NU) 

F( X5»XMU 

FOFX)  GO 

5)/2. 

-X4)/2. 

rXSIG 

IG/6. 

0 

X^)/2. 

-XM/2. 

EXSIG 

IG/6. 

0 

)  /2. 

(RNtNU) 
FCX5,XMU 
FOFX)  GO 
)/2. 

(RN,NU) 

F( X6,XMU 

FOFX)  00 

X4)/2. 

-X3)/2. 

FXSIG 

IG/6. 

0 

X3)  /2. 

-XD/2. 

CXSIG 

10/6. 

0 

X5 

(  P  M  ♦  N'  I ) 
F(X6,X*>U 
FOFX)  GO 


,XSIG) 
TO  9500 


,XSIG) 
TO  9250 


t  X  S I  0  ) 
TO  9125 


,XSIO) 
TO  9063 


, XSIO) 
TO  1313 


tXSIO) 
TO  1314 


,  X  S 1 0  ) 

TO  1316 


,XSIO) 
TO  1315 
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DELX=X5-X2 
FXMU=XB+0ELX/2. 
HXSIG=1.3*DELX 
GO    TO    7000 

1315    XB=Xt 

0ELX=X2-X6 
EXMU=XB-0FLX/4. 
cXSIG=**nELX 
GO  TH  700  0 

9125  X4=( X2  +  XI ) /2. 
N8R=NBR+1 
CALL  UNIF(RN,NU) 
FOFX=XNrDF(X^,XM 
IF(PM.LE.FOPX)  G 
X5=( X2+X3)/2. 
NBR=NBR+1 
CALL  UNIF(RN,NU) 
FOFX=XNCDF( X5, X^ 
IF(RN.LE.FOFX)  G 
X6=(B+X3) /2. 
NBR=MBR+1 
CALL  UNIF(PN,NU) 
FOFX=XN'CDF(  X6,X^1 
IF(RN.LF.FOFX)  G 
X7=2.*B-X6 

NBR=MBR+1 

CALL  UNIF(RN,NU) 

F0FX=XNCDF(X7tXM 

IFIPN'.LE.FO^X)  G 

XB=X7 

DFLX^XB-B 

EXMU=XP+0ELX/4. 

FXSIG=6.*DELX 

GO  TO  7000 
^012    XB  =  X3 

DELX=X6-XR 

EXMU=XB+DELX/2. 

FXSIG=1.3*0ELX 

GO  TO  7000 
9024    EXMU=( X3+X5) /2. 

EXSIG=(X3-Xr>)  '2. 

EXSIG=2.*FXS!G 

EXSIG=i=XSIG/6. 

GO  TO  7000 
9047    X6=( X4+X2) /2. 

NBR=NBR*-1 

CALL    UN  IF  (RN,  Nil) 

FOFX=XNCDF( X6,XM 

IFCPN.LE.FO^X)     G 

FXMU=( X5+X2) /2. 

FXSIG=(X5-X2) /?. 

EXSIG=2.*EXSIG 

EXSIG=FXSIG/6. 

GO  TO  7000 
9011  X7  =  2.*X4-X6 

N8R=NBR+1 

CALL    UNIF(RN»NUI 

F0FX=XNCDF(X7,XM 

IF(RM.LF.FOFX)     G 

XB  =  X6 

DELX=X?-XB 

FXMU=XB4-OELX/2. 

FXSIG=l.^nFLX 

GO  TO  7000 
9010    XB=X7 

DELX=X4-X7 

FXMU=XB-OELXM. 

EXSIG=6*0ELX 

GO  TO  7000 
9094    X5=2*X1-X4 

NBR=N8R+1 

CALL    UNTF(RM,NU) 


UtXSI 
0    TO 


UtXSI 
0    TO 


U  ,  X  S I 

0    TH 


U,XST 
0    TO 


G) 
909^ 


G) 

9047 


G) 

9024 


r,) 

9012 


U,XSI 
0    TO 


G) 
9011 


U  ,  X  S  T 

n   to 


G) 
9010 


30 


9003 


9250 


9006 


9005 


9004 


5555 


9007 


F0FX  = 
IF(RN 
XR  =  X4 
DELX  = 
EXMU  = 
EXSIG 
GO  IP 
XB=X5 
DELX  = 
EXMU= 
E  X  S  I G 
GO  TO 
X3=(£ 
NBR  =  N 
CALL 
FOFX  = 
IFCRN 
X4  =  (  X 
NBR  =  N 
CALL 
FOFX  = 
I F  (  R  M 
X5=(B 
NBR=M 
CALL 
FOFX  = 
IF(RN 
X6=2* 
NBR  =  N 
CALL 
FOFX  = 
IF(RN 
XB  =  X6 
DELX  = 
EXM|J= 
EXSIG 
GO  TO 
XB  =  X2 
DELX  = 
FXMIJ= 
EXSIG 
GO  TO 
EXMU= 
EXSIG 
FXSIG 
EXSIG 
GO  TO 
X5=(  X 
NBR  =  N 
CALL 
FOFX  = 
I F  (  R  NJ 
EXMU  = 
EXSIG 
EXSIG 
EXSIG 
GO  TO 
X6  =  2. 
NBR  =  N 
CALL 
FOFX  = 
IF(RN 
X8  =  X5 
OELX  = 
EXMIJ  = 
FXSIG 
GO  TO 
XB=X6 
DELX  = 
FXMU  = 
EXSIG 
GO    TO 


XNCOFCX 
.LE.FOF 

X2-X4 

XR+DFLX 
=1 .3*DF 
7000 

X1-X5 

XR-OELX 

=6.*DEL 

7000 
+X1)/?. 
BR  +  1 
UNIF(RN 
XNCOF(X 
.LE.FOc 
1 +X2)/2 
BR  +  1 
UMIF(RN 
XMCDF( X 
.LE.FOF 
+  X2) /2. 
RR  +  1 
UNIF(RN 
XNCDF( X 
.LE.FOF 
B-X5 
BR  +  1 
UNTFf&N 
XNCDF(X 
.LE.FOc 

X6-R 
XR+HELX 
=f*DELX 
7000 

X5-XR 

XB+DELX 

=1.3*DE 

7000 
( X2+X4) 
=( X2-X4 
=2.*EXS 
= EXSIG/ 

7000 
1+X3J/2 
RR  +  1 
UMIF(PNI 
XNCDF( X 
.LE.FOF 
( X1+X4) 
=(X4-X1 
=2.*FXS 
=  F  X  S  I  G  / 

7000 
*  X  3-  X  5 
RR  +  1 
l)NIF(PM 
XNCOF( X 
.LE.FOF 

Xl-X« 
XR+nFLX 
=1.3*DE 
7C00 

X3-XB 
XR-OFLX 
=  f>*OELX 
7000 


5,XMUt 
X)     GO 


/? 

LX 


X  S I  G  ) 
TH    QQ03 


/4. 
X 


»NU) 
3,XMUt 
X)     GO 


,NU1 
4,  XMI.I, 

X)     GO 


,  NU ) 
5,XV!Jt 
X)     GO 


6,XMU, 
X)     GH 


/4. 


/?. 

LX 

/2. 
)/2 

IG 
6. 


,NU) 

5»XM|Jt 

X)     GO 

/2. 

)/2. 

IG 

6. 


XSIG) 
TO    93  75 


XSIG) 

TO    9004 


XSIG) 
TO    c005 


XSIG) 
TO    QQ06 


XSIG) 
TO    5555 


,NUJ 

(S,  XMIJ. 
X)     GH 


/? 

LX 


XSIG) 
TO    9007 


M. 
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9375  X4=?.*A-X3 
NBR=NBR+1 

CALL  UNIF(PNfNU) 
FOFX  =  XNCDF<  X4,XMU» 
IF(PN.LE.FOPX)  GO 
XB=X3 

DELX=X1-XR 
FXMU=XR+OELX/2. 
EXSIG=1.3*OELX 
GO  TG  7000 

9376  XB=X4 
DELX=A-XB 
EXMU=XR-HELX/4. 
EXSIG=6*DELX 

GO  TO  7000 
9500  X2=(A+Xl)/2. 
NBR=NPR+l 
CALL  UNIF(RNfNU) 
F0FX=XNCDF(X2 tXMU, 
IF(RN.LE.FOFX)  GO 
X3=( Xl+B) /2. 
MBP=NRR+l 
CALL  UNIF(RNtNM) 
FOFX=XNC0F(X3,XVU, 
IF(RN.LE.FOFX)  GO 
X4=2.*R-X3 
NBR=NBR+1 

CALL  UNIF(RN,NU) 
FOFX  =  XNC0F{X4,XM|), 

IF(PN.LE.FOFX)  GO 
XB  =  X4 

DELX=X4-XR 
EX^U=XB+DELX/4. 
FXSIG=6*DELX 
GO  TO  7000 
5554  XR=Xl 

DELX=X3-XB 
EXMU=XB+DFLX/2. 
EXSIG=1.3*QELX 
GO  TO  7000 

5556  X4=( Xl+X2)/2. 
NBR=NBP+1 

CALL  UNIF(RN.NU) 

F0FX=XNCDF(X4tXM!}t 

IF(RN.LE.FOFX)  GO 

X5=( X1+X3) /2. 

NBR=NBR+1 

CALL  UNTF(RNfNU) 

FOFX=XNCDF( X5.XMU, 

IF(PN.LE.FT=X)  GO 

X6=2.*X3-X5 

NBR=NBR*l 

CALL  UNIF(RN9NU) 

FGFX=XNC0F(X6,XMUt 

IF(RN.LE.FOPX)  GO 

XB  =  X6 

DELX=XB-X3 

EXMU=XB+0ELX/4. 

EXSIG=6*DELX 

GO  TO  7000 

5559  XB=X1 

DELX=X5-Xl 
FXMU=XB+OELX/?. 
EXSIG=1.3*PELX 
GO  TO  7000 

5553  EXMU=( X4+X1) /2. 
EXSIG=(X]-X4)/?. 
EXSIG=FXSIG/3. 
GO  TO  7000 

5557  X5=(A+X2)/2. 
NBR=NBR+l 

CALL  UNIF(RN,NU) 


X  S I  G  ) 
TO  9376 


XSIG) 
TO  9501 


XSIG) 
TO  5556 


XSIG) 
TO  5554 


XSIG) 
TO  5557 


XSIG) 
TO  5553 


XSIG) 

TO  5559 
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FOFX=X 
IF(RN. 

X6=( X2 
NRR=NR 
CALL  ') 
FOFX=X 
I  F  (  R  N  . 
X7=(  X4 
NBR=NB 
CALL  U 
FOFX=X 
I F  (  P.  M  . 
XR=2.* 
NBR=NB 
CALL  U 
FHFX=X 
IF(PN. 
EXMU=X 
EXSIG= 
GO  TO 

123  FXMl)=( 
EXSIG= 
GO  m 

122  X«=(X6 
NBP  =  N'B 
CALL  U 
FO^X=X 
IF(RN. 
X9=( X^ 
NBR=NB 
CALL  U 
FOFX=X 
IF(RN. 
EXMU=( 
E  X  S I G= 
GO  TO 

125  FXMU=( 
EXSIG= 
GO  TO 

124  X9=(X2 
NBR=NB 
CALL  U 
FOFX=X 
I F  i  R  N  . 
EXMU=< 
EXSIG= 
GO    TO 

126  EXMU={ 
EXSIG= 
GO    TO 

121  X7=(X5 
NBR=NB 
CALL  U 
FOFX=X 
IF(RN. 
XR=( X6 
NBR=MB 
CALL  U 
FHFX=X 
I  F  (  R  N  . 
X9=( Xf 
NBR=NB 
CALL  U 
FOFX=X 
IF(RN. 
EXMU=( 
FXSIG= 
GO    TO 

129  EXMU=( 
FXSIG= 
GO    TO 

128       X9=(X7 


MCDF( X 

LE.FOF 

+  X4) fl 

R  +  l 

N I F  { R  N 

NCOF(X 

LF.FOf 

+  X1) /2 

R  +  l 

NIF(RN 

NCDr( X 

LE . F  OF 

X1-X7 

R  +  l 

NIFfRN 

NCDF(X 

LE.FOF 

«+( XR- 

6.*< X3 

7000 

X4+X7) 

1.3*(X 

7000 

+X4)/2 

R  +  l 

NIF1RNI 

NCOF(X 

LE.FOF 

+X7) /2 

R  +  l 

NIMRN 

NCDFCX 

LE.FOP 

X4+X9) 

1.3*(X 

7000 

X8+X4) 

( X4-XR 

7000 

+X6)/2 

R  +  l 

NIC(PN 

NCDF(X 

LE.FOF 

Xf+X8) 

( X8-X6 

7000 

X9+X6) 

l.?*(X 

7000 

+X2)/2 

R  +  l 

N  I  F  {  R  M 

NCDF(X 

LE.FOF 

+X2)/2 
R  +  l 

N I F  (  R  N 
NCDFCX 

LE.FOF 

+X4)/2 

R  +  l 

NTF( RN 

NCDF(X 

LE.FOF 

X6+X9) 

1 .3*( X 

7000 

XR+X6) 

( X6-X8 

7000 

+X2)/2 


,  X  S I G  ) 
TO    5561 


»XSIG) 
TO    121 


♦XSIG) 
TO    122 


,  X  S I G  ) 
TO    123 


5,  XMll 

X)     GO 


,NH) 
6,XMU 
X)    GO 


tMU) 
7fXM|.J 

X  )    GO 


R,XW|J 

X)  GO 

XI)  M 
-XI) 

/2. 

7-  X4 ) 


•  Nil) 

B.XMij,  XSTG) 
X)     GO    TO    124 


,NU) 
0  ,  XMM 
X)     GO 
/?. 
9-X4) 

/2. 
)/6. 


,NU) 
Q,  XMU 
X)     GO 
/2. 
)/6. 

r?. 

6-X9) 


»NU) 
7,  XWU 
X)     GO 


iNM) 
8  tXMij 
X)     GO 


,NH) 
9,  XMIJ 
X)     GO 
/2. 
9-X6) 

/2. 
)/6. 


»XSIG) 
TO    125 


tXSIG) 
TO    126 


,XSIG) 
TO    127 


,  X  S I G  ) 
TO    128 


,  X  S  I G  ) 
TO    129 
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CALL  UNTFCRN 
FOFX=XNCOF( X 
IFIRN.LF.FGP 
FXMU=(  X2+XP) 
EXSI0=( XR-X2 
GO    TO    7000 

130  EXMU=(  X9+X2) 
EXSIG=l.3*(X 
GO    TO    7000 

127  EXMH=( X7+X2) 
EXSIG=1.3*(X 
GO    TO    7000 

5561     X6=2.*A-X5 
NBR=NBR+1 
CALL    UMIFfPN 
FOFX=XNCDF(X 
IFIRN.LE.FOF 
XR=X5 
0ELX=X5-A 
EXMU=XR+OSlX 
EXSIG=1.3*05 
GO    TO    7000 

556C    XR  =  X5 

0ELX=A-X6 
EXMU=XB-OELX 
EXSIG=6*0tLX 
GO    TO    7000 

9501     X3=(A+X2)/2. 
NRR=NRR+1 
CALL    UN'IFfR.N 
FOFX=XNCDF(X 
IF(RN.LE.POP 
X4*< X1+X2) /? 
NBR=N&R+1 
CALL    UN  IF  CRN 
FOFX=XNCDF( X 
!F(RN.LE.FOF 
X5=2.*X1-X4 
N3R=NRR+1 
CALL     UNIF(RN 
FOFX=XNCDF( X 
IFCRN.LE.FOF 
XP  =  X^ 

DELX=XR-X1 
EXMU=XB+OELX 
EXSI0=6*DELX 
GO  TO  ^000 

9030  XR=X2 

DELX=X4-X2 
FXMU=XR+DFLX 
EXSIG-] .3*DE 
GO  TO  7000 

9504  X5=( X3+X2)/2 
NBR=NPR+1 
CALL  UNIF(RN 
FOFX  =  XN'COF(X 
IF(RN.LE.FO*= 
X6=( X2+X4) /2 
NBR=NBR+1 
CALL  HNIC(RN 
FOFX=XNCOF( X 
TF(RN.LF.FnP 

X7=2 .*X4-X6 

NBR=NRR+1 

CALL  UNIF(^i 

FOFX=XMCDF(X 

IF(RN.LE.FOF 

XP=X7 

DELX-X7-X4 

EXMU  =  XB*-PELX 

«rXSIG=6*DFLX 


i  Nil) 

9,  XMIJ, 

x)    on 

/2. 

)/6. 

/2. 
2-X°) 

/2. 

2-X7) 


,MII) 
6,X^Ut 
X)     GO 


XSIG) 

TO    130 


XSIG) 
TO    5560 


/2. 

LX 


/4. 


,NU) 
3, XMU, 
X)     GO 


tNU) 
4,XMMf 
X)     GO 


f  NU) 

5,XMU, 
X)    GO 


/4. 


/2. 

LX 


fMU) 
5,X*>U, 
X)     GO 


iNU) 

6.XMU, 

X)     GO 


XSIG) 
TO    9503 


XSIG) 
TO    9504 


XSIG) 
TO    9080 


tNUI 

7,X*MJ, 
X)    GO 


XSIG) 
TO    9081 


XSIG) 

TO    Q0R2 


XSIG) 
TO    9083 


/4. 
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9033 

9082 
9031 


9034 


^035 


9503 


9510 


9509 


9511 


GO    TO 

XB  =  X2 

D£LX=X 

FXM()=X 

EXSIG= 

GO    TQ 

EXMU=( 

FXSIG= 

GO    TO 

X6=( X3 

N'BR=N& 

CALL    U 

FOFX=X 

I  F  (  R  M  . 

FXMU=( 

F  X  S I  G= 

GO    Tn 

X7=2.* 

NBR=NR 

CALL    U 

FOFX=X 

IF(RN. 

X3=X6 

DELX=X 

FXMU=X 

FXSIG= 

GO    TO 

XB  =  X7 

DELX=A 

PXMU=X 

FXSIG= 

GO    TH 

X4=(X3 

NRR=NB 

CALL    U 

FOFX=X 

IF(PN. 

X5=< X2 

NBR=NB 

CALL    U 

FOFX=X 

IF(RN. 

X6=2.* 

NRR=NR 

CALL    U 

FOFX=X 

IF(PN. 

XR  =  X6 

DELX=X 

FXMU=X 

EXSIG= 

GO    TO 

XB  =  X3 

DELX=X 

FXMU=X 

FXSIG= 

GO    TO 

X6=( X3 

NRR=MB 

CALL    U 

FOFX=X 

IF(RN. 

FXMU=( 

FXSTG= 

EXSIG= 

EXSf 0= 

GO    TO 

EXMU=( 

FXSTn= 

FXSIG= 

EXSIG= 

GO    TO 


7000 

6-X2 

B  +  np 

1.3* 

7000 

X5  +  X 

(  X2- 

7000 

+  £  )  / 

R+l 

NIF( 

NCDF 

LE.C 

X3  +  X 

(  X*- 

7000 

A-X6 

R+l 

NIP{ 

MCD^ 

LE.F 


LX/2. 
DFLX 

2)/2. 

2. 

RMtNU) 
( Xfr.XWU 
HFX)  GO 
5)/2. 
X^)/6. 


,  X  S  T  G  ) 
TO  9034- 


RNtNU) 

(X7tXMU»XMG) 
OFX)  GO  TO  9035 


3-X6 

B+DFLX/2 
1.3*nELX 
7000 


-X7 
B-DF 
6*DE 
7000 
+  A)  / 
R  +  l 
NIF< 
NCDF 


LX/^. 
LX 

2. 


RNtNU) 

(X^,XMlJ 

LE.FOFX)  GO 
/2. 


+  X3) 

R  +  l 

NUF( 

NCOF 

Lf.F 

X2-X 

R  +  l 

NTF( 

NCOF 

LE.F 


RM,NIM 
(X5tXMU 
OPX)  GO 
5 

RN,NU) 
(X6.XMIJ 
OFX)  GO 


,  X  S  T  0  ) 
TO  9507 


, XSIG) 
TO  9509 


,  X  S I G  ) 

TO  9510 


6-X2 

B+OELX/6 
6*DELX 
7000 


5-X3 

B  +  DE 

1.3* 

7000 

+  X^) 

R  +  l 

NIF( 

NCDF 

LE.F 

X6+X 

(  X3- 

2.*c 

FXSI 

TOOO 

X4  +  X 

(  X6- 

2.*E 

CXSI 

7000 


LX/2. 
OELX 

/2. 

RN,NU) 
( X5,X^U 
OF  X )  GO 
3)  /2. 
X6)/2. 
XSIG 
G/6. 

6)  /2. 
X4)/2. 

XSIG 
G/5. 


.  X  S 1 0  ) 
TO  9511 
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9507 


9509 


7000 


1003 

1001 
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X5=( X4 
NRR=NR 
CALL  LI 
FOFX=X 
IF(RN. 
FXMU=( 
EXSIG= 
EXSIG= 
EXSIG= 

on  to 

EXMU=( 

FXSIG= 

EXSIG= 

EXSIG= 

GO  TO 

XR  =  0. 

OELX=0 

WRITF( 

IF(EXS 

OPT=0. 

DIF=0. 

XINC=0 

OPT=EX 

IF(OPT 

XTNC=( 

1NC=XT 

0IF=XI 

IF(DIF 

IF(DIF 

GO  TP 

A=OPT 

INC  =  0 

M=IO+I 

IS=70- 

M=IQ+I 


+  A) 
o  +  l 
NIP 
NCD 
LE. 
X^  + 
(  X4 
2.* 
FXS 
700 
X5  + 
(X4 
2.* 
FXS 
700 


/2. 

(RNtNU) 

F( X5,XMU, XSIG) 

FOcX)    GH    TG    9508 

X5)/2. 

-X5)/2. 

EXSIG 

IG/6. 

0 

A}/2. 

-X5)/2. 

FXSIG 

TG/6. 

0 


6, LOOM    FXMU»EXSIG 
IG.LT.O.)    EXSIG=-EXSIG 


MU-4.*EXSIG 

.LT.A)    GO    TG 

GPT-A)/FXSIG 

NC/i 

NC-INC 

.GE..5) 

.LT..5) 

1001 


NC  +  1 

NBR 
NC  +  1 


1003 


IMC=TNC+l 
TNC=INC 


CONDUCT    BRUCETON    TEST 


CLEAR    ARRAYS 

DO    10    1=1 t 200 
X(  I  )=0. 
IXO( I )=0 
IXX( I)=0 
NS( I )=0 
NG( I  1=0 
SUMAR-0. 
SUMRR=0. 
AR(I )=0. 
BR( I  )  =  0. 
10         CONTINUE 

LOAD    X    ARRAY 

DO    20    J=l,200 

X( J)=A+( J-1)*EXSIG 

20         CONTINUE 

CONDUCT    EXPERIMENT 
30         CMl    UNIP(RNtMU) 

FGFX  =  XNCDF(X(M)  »XMt)txSIG) 

IF(RN.GT.FGPX)GG  TO  40 

IXX(M)=IXX(M)+1 

M=M-1 

N=N+1 
5  0    IF(N.GT.IS)  GG  TG  60 

GO  TG  30 
40     IX0(M)=IX0(M)+1 

M  =  M+] 

N=N+1 

GO  TO  50 


PERFORM  BRUCETON  ANALYSIS 
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COUNT    RESPONSES    AMD    NON-RESPONSES 
6  0         ISUMX=0 
ISUMn=o 

no   i*  j=it2oo 

ISUMX=TSUMX*JXX( J) 
ISU^O=ISUM0+IXO( J) 
NS( J)=0 
AR( J )=0. 
RR( J  )  =  0. 
NG< J)=0 

14  CONTINUE 

DETERMINE  LESS  FREQUENT  EVENT  AND  LOAD  NS 
IF(  1SUMX.LE.ISUM0)  r,n  TO  15 
NT=!SUMO 
IFLAG-0 

00  21  J=l,200 
NS( J)=IXO( J) 

21  CONTINUE 
GO  TO  16 

15  NT= I  SUM X 
IFLAG=] 

DO  22  J=lt200 
NS( J)=TXX< J) 

22  CONTINUE 

DETERMINE  FIRST  AND  SECOND  MOMENTS 

16  JC0UNT=1 

17  IF(NS( JCOUNT).GT.O)  GO  TO  IS 
JC0UNT=JC0UNT*1 
!F(JCOUNT.GF.2O0)  GO  TO  104 
GO  TO  17 

IB    MC0UNT  =  200- .IC.OHNT 

DO  19  J  =  1  ,  MCOU'«'T 

NG( J)=NS( JCHUNT4 J-l ) 

AR( J)=( J-1)*NG( J) 

SUMAR=SUMAR+AR( J) 

BR( J  )=<<  J-l  )**2)*NGU) 

SUMPR  =  SUMRR4-B«  (  J) 
19    CONTINUE 

YPRIME=X( JCOUNT) 

CALCULATE    ESTIMATES    Oc    MEAN    AND    STANDARD    DEVIATION 

IE(  IFl  AG.FO.O)XMUEST  =  YPPIMF  +  EXSIG*(  (  SUMAR/NT  )  +  (  \m  /?. .  )  ) 
IE(.NOT.IFLAG.EO.O) XMUEST= YPR I M^+EXS IG* ( ( SUMAR/NT )-( 1 . 
SIGFAC=( (NT*Sl)MRR)-(SUMAR**2) ) /(NT**2) 
IF(SIGFAC.GT..3)     GO    Tn    1000 
EMU(LCOUNT)=0. 
EOEV(LCOUNT)=0. 
N0GP  =  N0G0«-1 
GO    TO    104 
1000    DEVEST=1.62*EXSIG*(SIGFAC+.029I 

LOAD    EMU    AND    EDEV 

FMU(LCRUNT)=XMUEST-XMU 
FDEV(LCOUNT)=DFVEST-XSIG 
ADDMU=ADDMU-»-FMU<  LCOUNT) 
ADDSIG=AnDSIG  +  EnEV( I  COUNT) 
ADDMUQ=ADDMUO+FMU( LC^UNT ) **? 
ADDSDO=ADDSDO+FDEV( LCOUNT) **2 
IF(PMU(LCOUMT)  .l.T.0.  )     GO    TO    91 
IF(EMU(LCniJNT)  .FO.O.)     GO    TO    92 
IMUHI  =  IMIJHI  +  1 

1  F  (  F Mll  (  L  COUM T  )  .  G T  . H  I  MM )     H I  Ml  j  =  F  MM  (  L C OUNT  ) 
IFl  .NO'.FMU(LCOUMT)  .GT.HIMU)     HI  MU  =  HI  mij 
GO    TO    93 

92         N0MU=N0MU4-1 

GO    TO    o^ 
91  IMUL0=IMUL0+1 

IF(EMU(LCOUNT)  .I.T.SMLO)     SMLO=EMU(LCOUNT) 
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)  DFVHI=FDFV(LCOUNT) 
OFVHI)  nEVHl=DEVHI 


IF(  .NOT.EMU(LCOUNT)  .  LT.SMl.O)  SMLO=SMLO 
^3    IF(FOEV(LCni)NT)  .IT.O.  )  GO  Tn  9^- 

IF(FDEV(LCDUNT) .EO.O.)  GO  TO  «5 

IDEVHI=IDEVHI+1 

IF(EDPV(LCnUMT) . GT.DEVHI 

IF(  .NOT.EDEV(LCOUNT)  .GT, 

GO  TO  104 
95    N0DEV=N0DEV+1 

GO  TO  )0<> 
94    IDEVL0=IDEVL0+1 

IF(LDEV(LCO'JNT).LT.DEVLO)  DEVLO=EDFV(LCOUNT ) 

IF( . NOT. FOE V(LC HUNT ).LT.DEVLO)  OEVLO=DFVLO 
104   JCOUNT=0 

SIGFAC=0. 

XMUFST=0. 

DEVEST=0. 

SUMAR=0. 

SUwBR=0. 

LC0UNT=LC0UNT+1 

HAVE    1000    EXPERIMENTS    BFEN    CONDUCTED    ? 

IF(LCOUNT.LT.lOOl)    GO    to    103 

TF  EXPERIMENTS  COMPLETED  CALCULATE  AND  WRITE  RESULTS 

EXN0GO=N0GO 

SAMAVM=ADDMU/( 1 000. -c XNOGO ) 
SAMAVD=ADDSIG/( 1000. -E X NOG H) 
SAMSOM=AODMUO/( Q 09. -F XNOGO) 
S AMSOD= ADD S^0/( 99°. -c XNOGO) 
IF( IANY.EO.l)  GO  TO  35 
IANY=IANY+1 
GO  TO  6«? 
35    STOP 
END 


SUBROUTINE  UNIF(RNtNU) 


SUBROUTINE  RETURNS  RANDOM  NUMBER  UNIFORM  ON  (0,1) 

REAL  MOD 

MOD=  2**31 

NR=129*NU+1 

RN=NR/MOD 

IF(RN.LT.O.O)  RN=-RN 

NU=NP 

RETURN 

END 


FUNCTION  XNCDF( V,XMU»SX) 


FUNCTION  SUBPRnr7pAM  CALCULATES 
X  IS  AN  R.V.  WITH  MEAN,XMijtAND 

ARG=( V-XMU)/SX 

XNCDF=SNCDF( ARG) 

RETURN 

END 


CUMULATIVE  NORMAL. 
STANDARD  DEVIATION, SX 


FUNCTION  SNCDF(X) 
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FUNCTION    SUBPROGRAM    CALCULATES    STANDARD    CUMIMATTVF    NORMAL. 
DATA    TPST/O.O/ 
IF( TEST. NE. O.O)    GO    TO    100 
SR2=    SQRT(P.O) 
TFST=1. 
100       SNCDF=( 1.0+ERF(X/SR2>) /2.0 
Rc  TURN 
END 
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